function [ExtFrac, aluminumOutIndex, cpxOutIndex,TransitionSp2Plag, TransitionGar2Sp, aluminumOutFrac, cpxOutFrac,MeltExtractedPb, ...
    BAI3, XOW_OUT, rXOW,initialBulkComp,FirstMeltPb,FirstAt2percent,FirstMeltIndex,FirstAt2Index,LastMeltPb,LastMeltIndex,...
    crustalthickness,poolCL_Iv,poolCL,poolBAI4_Iv,poolBAI4,running_mBAI4,running_CL_OUT,runningPressure] = ...
    PETROGEN2019_post_proc(Pb, T, ln, Lnexp, WnS, Wnr, Lnr, XM_OUT,BAI3, XOW_OUT, rXOW,initialBulkComp,CumFrac,...
    Gar2SpTransition_index, Sp2PlagTransition_index, Iv, CL_OUT, BAI4, rfracCrit,Dmatrix,CL_BAI3)


if rfracCrit == 1
    Lnexp = ln;
    BAI4 = BAI3;
    CL_OUT = CL_BAI3;
end


for kl = 1:size(Lnexp,2)
    Lnexp_cumulative(kl) = sum(Lnexp(1:kl)); % using (1-Wnr) has floating point issues, so use this instead
end
Lnexp_cumulative(Lnexp==0)=0;





ExtFrac = Lnexp_cumulative;
if size(Pb,2)>Gar2SpTransition_index
    Gar2SpTransition = Pb(Gar2SpTransition_index);
    Sp2PlagTransition = Pb(Sp2PlagTransition_index);
else
    Gar2SpTransition = NaN;
    Sp2PlagTransition = NaN;
end

%aluminumOutIndex = min([min(find(XM_OUT(:,6)==0 & Pb'>=Gar2SpTransition)) min(find(XM_OUT(:,5)==0 & Pb'<Gar2SpTransition & Pb' > Sp2PlagTransition))  min(find(XM_OUT(:,4)==0 & Pb'<Sp2PlagTransition))]);
cpxOutIndex = min(find(XM_OUT(:,1)==0));

PbRound = round(Pb,3);
aluminumOutIndex = min([min(find(XM_OUT(:,6)==0 & Pb'>Gar2SpTransition)) min(find(XM_OUT(:,5)==0 & Pb'<Gar2SpTransition & Pb' > Sp2PlagTransition))  min(find(XM_OUT(:,4)==0 & Pb'<Sp2PlagTransition))]);


aluminumOutFrac = CumFrac(aluminumOutIndex);
cpxOutFrac = CumFrac(cpxOutIndex);
MeltExtractedPb = Pb(find(ExtFrac >0,1));
FirstMeltPb = Pb(min(find(ln >0)));
FirstAt2percent = Pb(find(ExtFrac >0,1)); %Pb(min(find(CumFrac >0.02)));

FirstMeltIndex= (min(find(ln >0)));
FirstAt2Index = (min(find(ExtFrac >0,1)));



LastMeltPb = Pb(max(find(ln >0)));
LastMeltIndex = (max(find(ln >0)));


% if aluminumOutIndex < LastMeltIndex
%     LastMeltPb = Pb(aluminumOutIndex);
%     LastMeltIndex = aluminumOutIndex;
%     Lnexp(aluminumOutIndex:end) = 0;
% end


pooledliquid = (Lnexp*BAI3(:,1:9))./sum(Lnexp*BAI3(:,1:9)).*100;
if any(FirstAt2percent)>0
    crustalthicknessinkbars_novolumechange = FirstAt2percent.*max(ExtFrac); %if you change the mass from 1, divide by the inital mass
else
    crustalthicknessinkbars_novolumechange=0;
end

crustalthickness= crustalthicknessinkbars_novolumechange;


[difference index] = min(abs(Pb - crustalthicknessinkbars_novolumechange));


DEPTH=FirstAt2percent-LastMeltPb;
% %crustalthickness = FirstAt2percent.*initialdensity./averageliquiddensity(2)./(1+(1-max(ExtFrac))./max(ExtFrac));
% crustalthickness = DEPTH.*initialdensity./average_pooledliquid_density_novolumechangePT./(1+(1-max(ExtFrac))./max(ExtFrac));
% crustalthickness = crustalthickness./2;


%can probably comment out
% TransitionSp2Plag = ExtFrac(find(PbRound==Sp2PlagTransition));
% TransitionGar2Sp = ExtFrac(find(PbRound==Gar2SpTransition));



if size(Pb,2)>Gar2SpTransition_index
    TransitionSp2Plag = CumFrac(Sp2PlagTransition_index);
    TransitionGar2Sp = CumFrac(Gar2SpTransition_index);
else
    TransitionSp2Plag  = NaN;
    TransitionGar2Sp = NaN;
end



PhaseOut_ExtFrac = ExtFrac(aluminumOutIndex);


BAI3(:,10) =  ((BAI3(:,6)/40.3044))./(((BAI3(:,6)./40.3044)+((BAI3(:,5)./71.8464))));
XOW_OUT(:,10)= ((XOW_OUT(:,6)/40.3044))./(((XOW_OUT(:,6)./40.3044)+((XOW_OUT(:,5)./71.8464))));
rXOW(:,10) =((rXOW(:,6)/40.3044))./(((rXOW(:,6)./40.3044)+((rXOW(:,5)./71.8464))));

BAI3(:,11) =  BAI3(:,7)./BAI3(:,3);
XOW_OUT(:,11)= XOW_OUT(:,7)./XOW_OUT(:,3);
rXOW(:,11) = rXOW(:,7)./rXOW(:,3);

BAI3(:,12) =  (BAI3(:,8)+ BAI3(:,9))./(BAI3(:,8)+ BAI3(:,9)+BAI3(:,7));
XOW_OUT(:,12)= (XOW_OUT(:,8)+ XOW_OUT(:,9))./(XOW_OUT(:,8)+ XOW_OUT(:,9)+XOW_OUT(:,7));
rXOW(:,12) = (rXOW(:,8)+ rXOW(:,9))./(rXOW(:,8)+ rXOW(:,9)+rXOW(:,7));

initialBulkComp(:,10) = ((initialBulkComp(:,6)/40.3044))./(((initialBulkComp(:,6)./40.3044)+((initialBulkComp(:,5)./71.8464))));
initialBulkComp(:,11) = initialBulkComp(:,7)./initialBulkComp(:,3);
initialBulkComp(:,12) = (initialBulkComp(:,8)+ initialBulkComp(:,9))./(initialBulkComp(:,8)+ initialBulkComp(:,9)+initialBulkComp(:,7));


%% accumlating melt composition


if any(LastMeltPb) == 0
    LastMeltPb = 0;
end


M = Lnexp;
% Calculate melt productivity and crustal thickness
Mr = Iv.*(M');  % Melt production rate c(1/s)
Mr(Mr<0) = 0;      % Eliminate values where melt productivity is negative
Mr_norm =  Mr./sum(Mr);

poolCL_Iv = Mr_norm'*CL_OUT;        %weighted by melt velocity
poolCL = (Lnexp./sum(Lnexp))*CL_OUT; % NOTweighted by melt velocity

poolBAI4_Iv = Mr_norm'*BAI4;        %weighted by melt velocity
poolBAI4 = (Lnexp./sum(Lnexp))*BAI4; % NOTweighted by melt velocity


%[testindicies ,testindiciesb]= find(Lnexp);
[testindicies ,testindiciesb]= find(Mr');

for sss = 1:size(testindiciesb,2)
    Mr_here = Mr';
    running_mBAI4(sss,:)  = (Mr_here(1:testindiciesb(sss))./sum(Mr_here(1:testindiciesb(sss))))*BAI4(1:testindiciesb(sss),:);
    running_CL_OUT(sss,:)  = (Mr_here(1:testindiciesb(sss))./sum(Mr_here(1:testindiciesb(sss))))*CL_OUT(1:testindiciesb(sss),:);
    runningPressure(sss) = (Mr_here(1:testindiciesb(sss))./sum(Mr_here(1:testindiciesb(sss))))*Pb(1:testindiciesb(sss))';
    
end

if isempty(testindicies)==1
    running_mBAI4=[];
    running_CL_OUT=[];
    runningPressure=[];
end

if any(Lnexp) == 0
    running_CL_OUT = [];
    running_mBAI4 = [];
    BAI4_Mg8_running_shallow = [];
    BAI4_Mg8_running = [];
    BAI4_Mg8 = [];
    runningPressure = [];
end


end



